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> Abstract 

The ground-state energies of systems containing up to twelve 7r + 's in a spatial volume V ~ (2.5 fm) 3 
f — are computed in dynamical, mixed-action lattice QCD at a lattice spacing of ~ 0.125 fm for 

four different values of the light quark masses. Clean signals are seen for each ground state, 
allowing for a precise extraction of both the tt + tt + scattering length and 7r + 7r + 7r + -interaction from 
qq a correlated analysis of systems containing different numbers of 7r + 's. This extraction of the tt + tt + 

scattering length is consistent with than that from the 7r + 7r + -system alone. The large number 
J> of systems studied here significantly strengthens the arguments presented in our earlier work and 

unambiguously demonstrates the presence of a low energy 7r + 7r + 7r + -interaction. The equation 
of state of a ir + gas is investigated using our numerical results and the density dependence of 
the isospin chemical potential for these systems agrees well with the theoretical expectations of 
leading order chiral perturbation theory. The chemical potential is found to receive a substantial 
contribution from the 7r + 7r + 7r + -inter action at the lighter pion masses. An important technical 
aspect of this work is the demonstration of the necessity of performing propagator contractions in 
greater than double precision to extract the correct results. 
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I. INTRODUCTION 



Multi-hadron systems, from the deuteron to heavy nuclei to neutron stars, represent a 
significant fraction of the universe that we observe, and for decades the phenomeno logical 
study of these systems defined the field of nuclear physics. Understanding how nuclei and 
nuclear interactions emerge from Quantum Chromodynamics (QCD), the underlying theory 
of the strong interaction, is now a central goal of modern sub-atomic physics. Since hadrons 
are bound-states of quarks and gluons, they do not arise at any finite order in perturbation 
theory and a description from QCD has proved elusive. The only known non-perturbative 
method that systematically implements QCD from first principles is its formulation on a 
discretized space-time, lattice QCD. While it is still not possible to directly calculate the 
properties of even the simplest nucleus from QCD, the tools and technology are gradually 
being put in place to make such calculations possible with lattice QCD in the near future. 
The impact of the successful realization of this goal cannot be overstated. For the first time, 
it would allow reliable calculations of strongly-interacting many-body processes that are 
not (or are only poorly) accessible experimentally. Important examples are hyperon-nucleon 
interactions that many play significant roles in the interior of neutron stars. Further, it would 
enable the exploration of how the properties of such systems depend upon the fundamental 
constants of nature, exposing the (possible) fine-tunings between the light-quark masses that 
give rise to the multiple fine-tunings observed in nuclear physics. 

Our current understanding of nuclei requires a small but non-zero three-nucleon interac- 
tion. In the study of the structure of nuclei, we now have refined many-body techniques, such 
as Green's function Monte-Carlo (GFMC) [1J with which to calculate the ground states and 
excited states of light nuclei, with atomic number A < 14. Using modern nucleon-nucleon 
potentials that reproduce all scattering data below inelastic thresholds with x 2 /d°f ~ 1) 
such as AVi8 0, one fails, quite dramatically, to recover the structure of light nuclei. The 
inclusion of three-nucleon interactions greatly improves the predicted structure of nuclei, but 
at present, such interactions are difficult to constrain. At some point in the future, lattice 
QCD will be able to predict the interactions between multiple neutrons (and proton and 
mixed proton-neutron systems), bound or unbound in the same way it will be used to de- 
termine the two-body scattering parameters. A calculation of the three-neutron interaction, 
for instance, will be possible. 

As a first step toward the study of nuclei and their interactions using lattice QCD, in this 
work we study systems composed of up to twelve 7r + 's. These multi-pion systems are con- 
ceptually and computationally the simplest multi-hadron systems that can be constructed. 
In addition to the two-body interactions, there are expected contributions from multi-body 
interactions. Such multi-body interactions are not forbidden by the symmetries of QCD, and 
are expected to be present with a magnitude that can be estimated with naive dimensional 
analysis (NDA) [3J. However, they are qualitatively (and obviously quantitatively) different 
from systems involving nucleons. The lowest-lying continuum state of multiple pions in 
a large volume is perturbatively close to each pion carrying zero-momentum, whereas the 
wavefunctions of systems of multiple nucleons are subject to the Pauli exclusion principle. 

In a recent letter pQ, we have reported results of the first many-pion calculation 1 in lattice 
QCD. We explored systems containing up to five pions, extracted the n + n + scattering length, 



"Many-body" implies systems containing more than two bodies. When we refer to the pion we will mean 
the 7r + unless otherwise stated. 
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and, for the first time, found indications of a non-zero renormalization group invariant (RGI) 
7r + 7r + 7r + -interaction. Here we continue our investigations, presenting more detailed results 
of the previous studies and extending our work to systems containing up to twelve charged 
pions. We also use the recently derived expression for the ground-state energy of n-identical 
bosons in a finite volume at (D(L~ 7 ) [5] in our analysis, one order beyond that at which our 
previous calculations [I] were analyzed. 

Multi-pion systems are of interest in their own right as a strongly interacting boson gas at 
finite density and temperature. Such systems may be important for the late-time evolution 
of heavy ion collisions, such as those at RHIC and also in the interior of neutron stars (6J [7] . 
During the last several years there have been a number of theoretical explorations of pionic 
systems at finite isospin chemical potential, with or without a finite baryon number chemical 
potential. Leading order (LO) chiral perturbation theory (x?T) has been used to study the 
vacuum realignment that takes place in the presence of an isospin chemical potential that 
exceeds the mass of the pion, leading to a charged pion condensate, and excitations about 
this realigned vacuum. One of the important results of this present work is the calculation 
of the isospin chemical potential as a function of the isospin density. Our results are in good 
agreement with the LO %PT result, and further, demonstrate the sizable contribution from 
multi-pion interactions even at moderate densities. 

The structure of this paper is as follows. In Section [IT] we review the theoretical ex- 
pectations for the ground state energy of multi-pion systems at finite volume and discuss 
methods for extracting their interactions. In Section III we provide details of our lattice 
QCD measurements and analysis and in Sections [TV] and [V] we present the main results of 
our calculations. In Section |VI] we discuss the implications of our results for the equation 
of state of the pionic gas, its isospin chemical-potential and its pressure. Finally in Section 



VII we discuss the results in a global context and conclude. Certain technical details of 



contractions and numerical implementation are relegated to the Appendices. 



II. MULTI-MESON ENERGIES : ISOLATING THE TWO- AND THREE-BODY 
INTERACTIONS 

In recent works [51 El IS], the analytic volume dependence of the energy of n identical bosons in 
a periodic volume has been computed to 0(L~ 7 ), extending the classic results of Bogoliubov 
[TU] and Lee, Huang and Yang [TT]. The resulting shift in energy of n particles of mass M 
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due to their interactions is [5] 
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where the parameter a is related to the scattering length 2 , a, and the effective range, r, by 
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K. = 8.4019240, 
T x = 450.6392, (3) 



and n C m = n\/m\/(n — m)\. The three-body contribution to the energy-shift given in 

eq. ([l| is represented by the parameter ¥)% , which is a combination of the volume-dependent, 
renormalization group invariant quantity, r\\ , and contributions from the two-body scattering 
length and effective range, 
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The quantity 773 (//) is the coefficient of the three-7r + interaction that appears in the effective 
Hamiltonian density describing the system [5]. It is renormalization scale, /1, dependent. 
The quantity 5 is renormalization scheme dependent and we give its value in the minimal 
subtraction (MS) scheme, 5ms — —185.12506. 



2 In this work we use the Nuclear Physics sign convention for the scattering length, which is opposite to 
that of the Particle Physics sign convention. In this convention, the tt + tt + scattering length is positive. 
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For n = 2, the last two terms in eq. Q vanish and the remaining terms constitute the 
small a/L expansion of the exact eigenvalue equation derived by Lusher [121 [13] : 



pcot 5(p) = — S , (6) 




which is valid below the inelastic threshold. Below this threshold p cot 5(p) is the real part 
of the inverse scattering amplitude. The regulated three-dimensional sum is 



|j|<A 1 

s (*) = Eptp^ - 47rA ' (7) 

. |J| X 

where the summation is over all triplets of integers j such that | j | < A and the limit A — > oo 
is implicit. For n = 3, eq. ([!]) reproduces the shift of the ground state energy of three 
identical bosons that was recently calculated by Tan [9]. 

Naively, one might have expected to be able to determine the two-body effective range 
parameter, r, by calculating the energies of systems with different numbers of pions. How- 
ever, given the scattering parameter redefinitions of eq. ([2]), this is clearly not possible. 3 
Instead, the effective range will be extracted by calculating the energies of the excited states 
of, for example, the n = 2 pion system at finite volume, the lowest lying of which is pertur- 
batively close to the state in which the pions both carry (back-to-back) one unit of lattice 
momentum, 2tt/L. 

It is useful to form various combinations of the many-body energies in order to isolate 
or eliminate the various contributions from two-body or three-body interactions. Further, 
important checks can be made regarding the convergence of the large volume expansion, as 
more particles are added into the volume. In particular, the combinations involving systems 
with n, m and two bodies 
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vanish at order 0(L ) while 
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vanishes at order 0(L ) for any n, m. 



3 Writing — 1/pcot 6 = a = a + a 2 rp 2 /2 + and evaluating it at the shifted energy of two particles in the 
volume at LO in the volume expansion, gives a = a + 2ira 3 r/L 3 + .... 
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The three-body interaction, r) 3 , can be eliminated by forming combinations of the many- 
body energies, allowing for various determinations of a. One such combination is 
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which allows for the scattering length to be extracted at N 3 LO (omitting the last set of 
square brackets) and N 4 LO in the large volume expansion. 

Similarly, the three-body interaction can be isolated from combinations of the many-body 
energies. One such combination formed from the n-body and the two-body energies is 
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which allows for r] 3 to extracted at two orders in the expansion. This, of course, can be 
straightforwardly generalized to other combinations of energies that may or may not include 
AE 2 . 

The ground-state energies at finite volume, given in eq. ([!]), have been computed in non- 
relativistic quantum mechanics, with relativistic effects added perturbatively j5]. They are 
given in terms of the scattering parameters, and the three-body interaction. For pionic 
systems in particular, it is natural to ask about the role of chiral perturbation theory in 
such a calculation. In xPT, the expansion parameters for this system (in the p-regime) 
are p/A x and m 7T /A x , where A x ~ 4:Tcf n is the chiral symmetry breaking scale. The 
expression for the ground state energy of n-7r's in a finite volume (which remain to be 
determined) will be a dual expansion in 27r/(A x L) and m n /A x . As such, in order to extract 
the three-body interaction that first enters at 0(L~ 6 ), the calculation in x?T would need to 
be performed at N 3 LO. While providing the complete quark-mass dependence at this order, 
such a calculation would involve the evaluation of three-loop diagrams, and contributions 
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from the relevant N 3 LO counterterms. Organizing the perturbative expansion using non- 
relativistic effective field theory, EFT(^) (see e.g. Refs. [151 US]); greatly simplifies the 
result, but does so at the expense of ignorance of the quark mass-dependence of the EFT(^) 
parameters without further calculation (such as the quark mass dependence of the scattering 
length) . 

It is important to notice that for the systems containing a large enough number of pions, 
the energy-shift of the ground state exceeds that required to pair-produce pions. The leading 
relativistic effects that are included in the analytic expression for the energy-shift include 
only the single particle relativistic kinematics and corrections to the two-pion scattering 
amplitude, they do not include contributions due to inelasticities, including pair-production. 
However, we conclude that such effects make a small contribution to the ground-state energy 
as we find no evidence for deviations from eq. ([I]). Nonetheless, this aspect of these systems 
must be explored further. 

III. METHODOLOGY AND DETAILS OF THE LATTICE CALCULATION 

A. Lattice configurations and quark propagators 

The results of the numerical computations presented in this paper were obtained using the 
mixed-action lattice QCD scheme developed by LHPC (13, HE] and are based on the coarse 
MILC lattice configurations [19]. These lattices have a lattice spacing of b ~ 0.125 fm, and 
a spatial extent of L ~ 2.5 fm. They were generated using the asqtad-improved [201 EI] 
staggered formulation of lattice fermions, taking the fourth root of the fermion determinant, 
and the one-loop, tadpole-improved Symanzik gauge action [22]. Herein, we assume that 
the "fourth-root trick" 4 recovers the correct continuum limit of QCD. These ensembles of 
configurations have a fixed (almost physical) strange quark mass while the degenerate light 
quarks were varied over a range of masses; see Table [J and Refs. [3S1 ESI SHI SB S2] for 
details. 

Based on these configurations, valence quark propagators using the domain-wall (DW) 
formulation of the lattice fermion action (131 SH SHI SEl SZ] were computed from smeared 
sources on each gauge-field configuration. Hyper-cubic (HYP) smearing [2H SHI SSI 150] 
was applied to the gauge links used in the domain-wall fermion action to improve chiral 
symmetry, and in calculating the quark propagators, Dirichlet boundary conditions were 
imposed to reduce the original temporal extent of 64 down to 32. This procedure is optimized 
for nucleon physics and indeed leads to minimal degradation of a nucleon signal, however it 
does limit the number of time slices available for fitting meson properties in which the ratio 
of signal to noise remains constant in time. Further details about the mixed-action scheme 
can be found in Refs. [381 El]- A summary of the lattice parameters and resources used 
in this work is given in Table [Tj In order to generate large statistics on the existing MILC 
configurations, multiple propagators from sources displaced both temporally and spatially 
on the lattice were computed. 

In the continuum chiral limit the rif = 2 staggered action has an SU(8) l®SU (8)#Cg>?7(l)y 
chiral symmetry due to the four-fold taste degeneracy of each flavor, and each pion has 15 



4 For an introduction to staggered fermions and the fourth-root trick, see Ref. [23 . For the most recent 
discussions of the topic, see Ref. [M1I13I2S1I1EI1H1[M1[^[3I1[M1IM1IM1IM1IM1[S7] 



7 



TABLE I: The parameters of the MILC gauge configurations and domain-wall propagators used in 
this work. The subscript I denotes light quark (up and down), and s denotes the strange quark. The 
superscript dwf denotes the bare-quark mass for the domain-wall fermion propagator calculation. 
The last column is the number of configurations times the number of sources per configuration. 
Throughout the paper we will use the pion mass to refer to the ensembles. 



Ensemble 


bmi 


bm s 


l dw f 

bm^ 


j dw f 

bm s 


[ MeV ] 


# of propagators 


2064f21b676m007m050 


0.007 


0.050 


0.0081 


0.081 


~ 291 


1039 


X 


24 


2064f21b676m010m050 


0.010 


0.050 


0.0138 


0.081 


~ 352 


769 


X 


24 


2064f21b679m020m050 


0.020 


0.050 


0.0313 


0.081 


~ 491 


486 


X 


24 


2064f21b681m030m050 


0.030 


0.050 


0.0478 


0.081 


~ 591 


564 


X 


23 



TABLE II: The mass and decay constant of the ir + , and the energy-shift and scattering parameters 
in the tt + tt + system calculated previously |54J . The first uncertainties are statistical, the second 
uncertainties are systematic uncertainties due to fitting and the third uncertainty, when present, is a 
comprehensive systematic uncertainty [53]. Cvr is the one-loop counterterm in contributing 
to 7r + 7r + scattering, and 5 is the phase-shift. Recall that we are using the nuclear physics convention 
for the sign of the scattering length. 



Quantity 


m„ ~ 291 MeV 


m w ~ 352 MeV 


m„ ~ 491 MeV 


m-n ~ 591 MeV 


Fit Range 


8-12 


8-13 


7-13 


9-12 


Tfl-it (l.U.) 

/» (Lu.) 
rrhrl h 


0.18454(58) (51) 
0.09273(29) (42) 
1.990(11)(14) 


0.22294(31) (09) 
0.09597(16)(10) 
2.3230(57)(30) 


0.31132(28)(21) 
0.10179(12)(28) 
3.0585(49)(95) 


0.37407(49)(12) 
0.10759(28)(17) 
3.4758(98) (60) 


Fit Range 


11 - 15 


9-15 


10 - 15 


12 - 17 


AE 2 (Lu.) 
m^a^T (6 / 0) 

lii =2) {b + 0) 
S (b ^ 0) (degrees) 
\p\/m n 


0.00779(47) (14) 
0.1458(78) (25) (14) 
6.1(1.9)(0.7)(0.4) 
-1.71(14)(04) 
0.2032(60)(18) 


0.00745(20)(07) 
0.2061(49)(17)(20) 
5.23(68)(24)(28) 
-2.181(81)(28) 
0.1836(25)(09) 


0.00678(18)(20) 
0.3540(68)(89)(35) 
6.53(32) (42) (16) 
-3.01(09)(12) 
0.1480(17) (23) 


0.00627(23)(10) 
0.465(14) (06) (05) 

6.90(40) (18) (13) 
-3.46(17)(07) 
0.1298(24) (10) 



degenerate partners. At finite lattice spacing, this symmetry is broken and the taste multi- 
plets are no longer degenerate, but have splittings that are (D(a 2 b 2 ) for the asqtad staggered 
action. When determining the mass of the DW valence quarks there is an ambiguity due 
to the non-degeneracy of the 16 staggered bosons associated with each pion. One could 
choose to match to the taste-singlet meson or to any of the mesons that become degenerate 
in the continuum limit. The choice of tuning to the lightest taste of staggered meson mass, 
as opposed to one of the other tastes, provides for the "most crural" domain-wall mesons 
and therefore reduces the uncertainty in extrapolating to the physical point. The mass 
splitting between the domain-wall mesons and the staggered taste-identity mesons, which 
characterizes the unitarity violations present in the calculation, is then given by [52J [33] 

b 2 m 2 Wi - b 2 ml dwf = 0.0769(22) . (12) 

Simple properties of the tt + have been computed to high precision on the ensembles of 
coarse MILC lattices that are used in this work, both with staggered valence quarks and 
domain-wall valence quarks. Further, the energies of the lowest-lying 2-rr + states, which 
lead directly to the scattering lengths using Liischer's method, have been computed rela- 
tively precisely [54J. The results of the previous mixed-action calculations on these lattice 



8 



ensembles [51] are shown in Table [Tip . 

The results of the present calculation are presented in lattice units (l.u.), or in terms of 
dimensionless quantities such as m^/ f w which eliminates the requirement of scale setting. 
They are performed only at one lattice spacing, due to limited computer time, and as a 
result the continuum limit cannot be determined. Unlike the two meson system, for which 
mixed-action chiral perturbation theory (MA^PT) j55j [56J, [57] has been used to include the 
leading order effects of the finite lattice spacing, MA^PT calculations have not yet been 
performed for the multi-7r + systems, and therefore the leading lattice spacing artifacts in 
these calculations cannot be removed at present. The lattice spacing artifacts are assumed 
to be small, occurring at 0(b 2 ), but a systematic study must be performed in the future. 



B. Correlation functions 

In this work we determine the tt + tt + and 7r + 7r + 7r + interactions from the ground-state energy 
of n < 13 7r + 's (isospin stretched states). By working in the m u = rrid limit and restricting 
the calculation to states of maximal isospin, only the simplest sets of propagator contractions 
are required to be performed (i.e. no disconnected diagrams) in order to form the correlation 
functions from which the ground-state energies are extracted. 

Naively, there are (n!) 2 contractions (for large n this behaves as ~ 
contributing to the correlation function of n-7r + 's, 

C n (t) oc (^£>-(x,f)j (V(O,0)j ) , (13) 

where 7r + (x, t) = u(x, t)^d(x., t). However, this correlation function can be written as 6 

C n (t) oc ( ( rjUV T ) , (14) 

where 

n = S(x,t;0,0) 5 f (x,t;0,0) , (15) 

X 

and S(x, t; 0, 0) is a light-quark propagator. The object (block) II is a 12 x 12 (4-spin and 
3 color) bosonic time-dependent matrix, and rj a is a twelve component Grassmann variable. 
Using 

{rr^...rf n W,-rip n ) oc e***-*.*-**- e Pl p M . 4l2 _ n , (16) 
leads to correlation functions 

c n (t) = e^-^-^ e MMM2 _ n (n)£ (n)£ .. (n)£ . (17) 



5 Until this point the two-body scattering length for a generic system has been denoted by a. For the tt + tt + 

—(1—2) 

system, we denote the scattering length by . 

6 We thank David Kaplan and Michael Endres for discussions on this topic. For a general approach to 
evaluating contractions involving a large number of fcrmions, see Ref. [58 . 
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(c) 



o q> 

(d) 



FIG. 1: Graphical representation of the contractions for three pions with I z = 3, (a,b,c). By 
restricting to the maximal isospin, computationally demanding contractions such as the type shown 
in (d) are eliminated. 



While correct, further simplifications are possible. Let us recall that for an arbitrary 12 x 12 
<^ « l A.4 i — e PM2 (1 + \Af a \ (1 + \A)% . . . (1 + XAtl 



1 

12! 
1 

12! 



Ct\OL2-.OL\2 



\ L2.fi OL\OL2-.a\2 , 



a r 1 + 



+ A" 12 C n £^2..a n ^..ix2-n £/3i 



lh~Pn(,\~£.V2-n 



+ X £ 



12 aia2..ai2 



£ /3i/3 2 ../3i2 



.4 



.4 



A Y 1 

\PL2 
I a\2 



A 



A 



± J2 n C 3 A J Cj(t) 



(18) 



i=L 



where in the last line we identify the matrix A with n. Further, 



det (1 + XA) = exp (Tr [log [ 1 + XA] ] 

= 1 + ATr[A] + ^ 

A 3 



exp Tr 



E 

.p=i 



P 



-X P A P 



({Tr[A}) 2 - Tr [ A 2 ] ) 



+ — ( 2Tr [ A 3 ] - 3Tr [ A } Tr [ A 2 ] + ( Tr [ A ]) 3 ) + 



.(19) 



Therefore, by equating terms of the same order in the expansion parameter A in eq. (18) 



and eq. (19), one can recover the n-ir + correlation functions in eq. (17). As an example, the 



contractions for the 3-ir + system are 

C 3 (t) oc tr c ,s [IT] 3 - 3 tr 0lS [n 2 ] tr c ,s M + 2 tr c> s [n 3 ] 



(20) 



where the traces, trc,s, are over color and spin indices. The three contributions in the 
correlator in eq. (20) are shown in fig. [l], (a), (b), and (c), respectively. As it is the energy of 
states with maximal z-component of isospin that are calculated in this work, disconnected 
contractions, such as those in fig.fTld), do not contribute to the correlation functions that are 
computed. The explicit form of the contractions for n = 1, ... 13 are given in Appendix [A| 
Rewriting the contractions in terms of traces over the Il-blocks greatly reduces the required 
number of calculations, with the number of independent contributions to the correlation 
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function equal to the partition, P(n), of n objects. An estimate of the number of operations 
that must be performed to generate the correlator for n-mesons is ~ n(12 x 13 — 1) + 

YTj=i P(J)> which f° r large n scales as ~ ^/^~^ e7T ^^ using a classic result of Hardy and 
Ramanujan [59]. While for n = 12 there are ~ 2.3 x 10 17 independent contractions that must 
be performed, this can be accomplished with ~ 2 x 10 3 calculations to produce the ~ 80 
terms contributing to the contraction. Since, in this work, each contraction is performed 
with only a single quark propagator on each configuration, the Pauli-exclusion principle 
requires that the n > 13 identical meson contractions vanish identically, e.g. Ci3(t) — V t, 



implying the A 13 and higher terms in the expansion of eq. (19) vanish. Written in terms of 
contractions of propagators in flavor and color space, the n = 13 case of eq. (19), represents 
a generalized Cayley-Hamilton identity satisfied by all matrices of size less than 13 x 13. To 
perform calculations on systems containing more than twelve pions, additional propagators 
will be required. 

C. High-precision implementation 
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FIG. 2: The n = 13 correlation function as a function of the precision used to perform the 
calculation. The vertical axis is the logarithm of the sum over time-slices of the absolute value of 
the n = 13 correlation function on a representative gauge field configuration. 

In order to calculate the n-n + correlation functions, particularly for n > 8, it is necessary to 
use a numerical representation with precision greater than that of standard 64-bit machine 
precision. 7 This need arises because of the large products of propagators that must be 
computed and is not particular to the contractions studied here. In particular, the numerical 
issues impacting calculations of multi-pion systems that we have found in this work will also 
impact calculations of multi-nucleon systems. 



' The failure of double precision operations for these correlation functions is explored in detail in Ap- 
pendix |B| 
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Our implementation of these contractions uses arbitrary precision arithmetic based on 
the ARPREC library |60j which was extended for the particular operations needed here, ma- 
trix multiplications and traces. For the correlators studied here 64 decimal digit precision 
(approximately octupule precision) in internal operations is sufficient to give results accurate 
to sixteen digits. The additional overhead of using this numerical representation causes the 
high precision contraction code to run ~10-50 times slower than a double-precision version 
but is only marginally dependent on the precision used 8 . For the n — 13 correlation function 
it is instructive to look at the dependence of the resulting correlator on the precision used 
in the computations. Since the correlation function must vanish identically for any input 
propagator, it is a very stringent test of the codes used herein. In fig. |2j the logarithm of 
the sum over time-slices of the absolute value of the correlation function as a function of 
the digits of precision used to perform the contractions on a representative configuration is 
shown. From extrapolating the results shown here, we conclude that the correlator is indeed 
identically zero. 



D. Analysis 

The correlation functions from which we extract the ground-state energy of the n-7r + system 



are given in eq. (13), and, on a lattice with infinite extent in the time direction, behave as 

C n (t) ^ X n) e~ E - * (21) 

at large times. It is the difference between this energy, E n and n times the tt + rest mass 
that is equated to the energy difference given in eq. ([!]), and which is extracted from the 
ratio of correlation functions 

G n (t) = [ ||^ f ^ e -**> t , (22) 

where AE n is that of eq. ([!]). While there are a number of ways to extract the energy dif- 
ference from the correlation function, perhaps the most visually pleasing one is to construct 
the effective energy difference function, defined to be 



AK <" = to HG^Ti)J ™ AE " ' (23) 

In the limit of an infinite number of measurements, this function would tend to a constant 
equal to the ground-state energy splitting. Of course, for any real calculation, both the 
number of gauge fields and the number of propagators per gauge field are finite, and as such 
the object AE^f'(t) consists of a central value and an associated uncertainty at each time- 
slice, t. Further, the temporal extent of the lattice is finite, giving rise to both forward and 
backward propagating hadrons. As such, AE^- (t) is constant (up to statistical fluctuations) 
only over a finite number of time-slices, in a region between where the forward propagating 
excited states have died-out sufficiently, and where the backward propagating states have 



Checks of our C++ contractions have also been performed using Mathematica 6.0 (time costs prevent us 
using this on a large scale). 
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not yet become significant. In the current situation where Dirichlet boundary conditions are 
imposed, the behavior of the correlator is modified by the "reflection" of forward and back- 
ward propagating states from the boundaries. As these reflections are poorly understood, 
the region close to the boundary is omitted in our analysis. 

In our lattice calculations, multiple propagators and correlation functions are computed 
on each gauge configuration. These propagators and correlation functions are not statisti- 
cally independent unless the sources are separated by many correlation lengths. To account 
for this, all the correlation functions for a fixed n computed on a given configuration are 
averaged (blocked) into one correlation function. The configurations are found to be statis- 
tically independent 9 , and the blocked correlators on each configuration form the basis of our 
statistical analysis. 

In order to extract the energy difference AE n from AE^-(t), a fitting interval must be 
selected. This interval is chosen to be entirely contained in the region where the AEf-(t) 
is consistent with a constant. Once the fitting interval has been selected a correlated \ 2 
minimization is performed to extract the parameter AE n , defined in eq. (22). The covariance 
matrix that determines the correlated weightings of each of the values of AE^'it) on any 
given time-slice is generated using a single-subtraction Jackknife procedure. 10 The central 
value of AE n is the value that minimizes the correlated x 2 > an d the standard statistical 
uncertainty is determined by the values of AE n for which \ 2 ~^ X 2 + l- The fitting systematic 
uncertainty associated with the fitting procedure is determined by varying each end of the 
fitting range by —2 < At < +2, and refitting the energy-splitting. 

To study the scattering length and three-body parameter, rj 3 , using eqs. (10) and (11), 
the appropriate ratios of the C n {t) correlators are used to define effective scattering length 
functions for each n (LO, NLO, N 2 LO) or each pair {n,m} (N 3 LO and N 4 LO) and effective 
rj 3 functions for each n. We then analyze these in the same manner as the energy differences 
above. This leads to multiple determinations of a^ 2 ^ and r} 3 . The effective functions defined 
by C6,7j ec L s - @ an d @, are studied similarly. 

To make use of the full data set, we also perform a simultaneous, corre- 
lated fit of a^j" 2 ) and 773 to the effective masses for n = 2,3,...,N max for 
N ma , x = 3, . . . , 12. In order to do this, fitting ranges, tffi n < t < tmax, are 
chosen for each n (as above) and the data is assembled into a vector V = 
{A£f -(41), • • • AEfitgl), AEfit® n ), ■ ■ ■ AEfitgL), . . . AEfHt^), . . . AEg($&)}. 
A correlated \ 2 minimization is performed to extract the parameters a^n an d V3 v i a the 
fit vector U = {AE 2 , AE 2 , AE 3 , AE 3 , AE 12 , AE 12 }, where AE n is given in 
terms of a^ 2 ^ and r/3 using eq. (fxj) , and where the covariance matrix that determines the 
correlated weightings of each contribution is generated using the Jackknife procedure. The 
standard (statistical) uncertainties on a^" 2 ^ and fj 3 are determined from the maximum and 
minimum values of each parameter dictated by the uncertainty-ellipse corresponding to 
their values for which x 2 ~ > X 2 + 1- The systematic uncertainty associated with the fitting 
procedure is determined by repeatedly and randomly varying each end of the fitting range 



9 This is tested by averaging over sets of neighboring configurations and performing analysis on the resulting 
blocked ensemble. For block sizes of 1, 4 and 12, no noticeable difference is seen. 

As a check, we also performed a separate analysis using bootstrap re-sampling. The resulting energies 
and parameters were consistent, and, for simplicity, we focus on a single analysis in the main discussion. 
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for each correlator by —2 < At < +2, refitting the parameters a^ 2 -* and fj 3 , determining 
the complete range of values for each parameter associated with \ 2 ~^ X 2 + 1- The 
systematic and statistical uncertainties are combined in quadrature in this work. 



IV. LATTICE QCD RESULTS 

Using the techniques discussed in the previous section, we now turn to analysis of the results 
of the lattice calculations. Our main aim is to extract the parameters a^= 2 ^ and Tj 3 which 
we can do in a number of ways, either by forming particular combinations of energies, 



eqs. (10) and (11), or by a coupled analysis accounting for correlations among different n. 
The different methods give consistent results but we find that the most precise extraction 
is achieved using the latter method and consequently our final results are generated from 
this technique. Before we present these results, we first detail the simpler analysis using 
combinations of energies. As there are a large number of correlation functions that we study 
in this work, in some intermediate stages, we only display results for a single quark mass 
corresponding to m n = 291 MeV (in terms of uncertainties, this ensemble is neither the best 
nor the worst). 



A. Multi-pion energies and energy differences 

A priori, it may seem surprising that the correlation function of twelve 7r + 's can be calculated 
at all. On the ensemble associated with the lightest pion mass, m n = 291 MeV, the 12-7r + 
state has an energy of E 12 ~ 3.5 GeV, while on the ensemble with m n = 591 MeV, the 
12-7r + state has an energy of Eyi ~ 7.1 GeV. It is not immediately obvious that such a 
rapidly diminishing exponential can be cleanly measured but in these systems we find that it 
is possible. 12 As an example, for m n = 291 MeV, the effective energies (in lattice units) for 
systems with n = 1, 2, .., 12 7r + 's are shown in fig. [3j Well-defined plateaus in the effective 
energy plots are seen for all systems, with the relative statistical uncertainty in the data 
almost constant as a function of n. As can be seen from the fits to the energies (statistical 
and systematic uncertainties are shown in quadrature), the precision with which the energy 
can be extracted is high, typically < 2%. The total relative uncertainties on the energy of the 
n-Ti + systems are shown for all data sets in fig.[|J with only a slight dependence of n apparent. 
An additional point is that it was not obviously the case that the Gaussian-smeared source 
for the light-quark, that is suitable for the single pion ground-state, would have sufficient 
overlap onto the multi-pion ground-state to produce useful correlation functions. However, 
it is clear that it did. 

The energy differences which enter into eq. ([T| and subsequent results can also be ex- 
tracted cleanly, although with somewhat less precision than the individual energies. The 
effective energy difference plots, along with our fits to the energy differences are shown in 
fig. [5] (again for the = 291 MeV ensemble) while the relative uncertainties in our extrac- 



11 An alternative systematic procedure of repeatedly and randomly choosing triplets of time-slices in each 
fit range ±1 time-slice and refitting is also used, giving qualitatively similar results. 

12 In purely pionic systems there is no exponential degradation of the signal-to-noise ratio, unlike in most 
other hadronic systems [51] . 
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FIG. 3: Effective energy plots for the n = 1, . . . , 12 pion correlations as a function of time for the 
ensemble with = 291 MeV. The solid line and shaded region show the fitted energy and the 
systematic and statistical uncertainties combined in quadrature. 



tions are given in fig. |4} All effective energy splitting plots show behavior that is consistent 
with a single exponential (within statistical uncertainties) for a number of time slices. As 
discussed above, the region above t/b ~ 16 is contaminated by reflections from the Dirichlet 
boundary at t/b = 22 and is discarded in our analysis. 



B. 7r + 7r + scattering length 

Since the 7r + 7r + system is free from three- and higher-body hadronic interactions, it is the 
ideal place to extract the two body parameter, a^ 2 \ As is well known, this can be done 
without resorting to an expansion in a^^/L using the eigenvalue equation in eq. 
We refer to a^ 2 -* extracted in this way as the Liischer result and it forms a benchmark for 
extractions in the n > 2 systems. Eq. tun also allows us to extract a^ 2 ^ in a number of ways. 
At orders L~ 3 , L~ 4 , L~ 5 (LO, NLO and N 2 LO, respectively), each e nerg y difference, AE n for 



n 



2, . . . , 12, leads to a separate extraction of 2 >. Finally eq. (10), allows us to extend 



these extractions to N 3 LO and N 4 LO in the 2 - ) / L expansion by combining the n- and m 
body energy differences to eliminate three-body interactions. Choosing 3 < m < n < 12, 
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FIG. 4: Relative uncertainties in extractions of energies (left) and energy differences (right) as a 
function of the number of 7r + s. Results are shown for all quark masses. 



allows for 45 separate extractions. 

As a representative example, fig. [H] shows the LO, NLO and N 2 LO effective «i^r 2 ^ plots 13 
for n = 8, and the N 3 LO effective a^ 2 " 1 plots for {n, m} = {8,4} from eq. (JIOJ) , for the 
m n = 291 MeV ensemble. A summary of all the extractions of a^ 2 -* is given in fig. [7] (the 
N 4 LO extractions are entirely consistent with those at N 3 LO in all cases and are omitted). 
For the n = 2 data, it is clear that NLO and higher extractions, fig. [7| yields the same 
a^ 2 -* as the exact eigenvalue method of Liischer. However, for the multi-pion systems, an 
n-dependent systematic deviation from the exact eigenvalue method of Liischer is found at 
LO, NLO and N 2 LO. This is particularly clear for the lighter mass ensembles. In contrast, 
the extractions of a^ 2 ' in which the 7r + 7r + 7r + -interaction (or at least, a term that behaves as 
n Cz) is eliminated (N 3 LO and N 4 LO) are in close agreement with the n = 2 exact eigenvalue 
result for all n. From this alone we conclude that the calculation suggests the presence of a 
significant 7r + 7r + 7r + -interaction. 

We can further demonstrate the need for a 7r + 7r + 7r + -interaction by using eq. (jTl) to com- 
pute the energy shifts at 0(L~ 7 ) using the value of a^ 2 -* extracted from the 7r + 7r + -system 
using the exact eigenvalue method and setting fj 3 = 0. These can then be compared to the 
calculated effective energy differences as shown in fig. [8] for the m n = 291 MeV ensemble. 
The deviations between the predictions and the effective energy splittings are significant 
and grow with increasing n. Over the same sets of time-slices as used in the fits to the 
two-particle energy difference, the x 2 /d-o.f. of such an ansatz is 8.62 and therefore r] 3 = 
very poorly describes the results of the calculation. 



13 Given the accuracy with which m, and /„• have been extracted on these ensembles, the uncertainties in 
TnirCinT^ 2 ^ and m^f^rj^ are dominated by the uncertainties in aiir -2 ^ and rj^ , respectively. 
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FIG. 5: Effective energy difference plots for n = 2, ... ,12 pion correlations as a function of time 
for the ensemble at = 291 MeV. The solid line and shaded region show our fitted energy and 
the systematic and statistical uncertainties combined in quadrature. 



C. Three body interaction 



The 7r + 7r + 7r + -interaction can be explicitly constructed using eq. (11) for n — 3, . . . , 12. As 
an example, the effective fj^ plots for the m n = 291 MeV ensemble are shown in fig. [9J A 
clear plateau inconsistent with zero is seen in most cases. Fig. 10 shows the n dependence 
of the extracted value of 7/3 for each quark mass. In general, the combined systematic and 
statistical uncertainty of the extractions decreases with increasing number of 7r +, s. This is 
not surprising given the combinatoric factors that appear in the expression for the energy 
shift. 



D. Convergence: ^6,7 

Before presenting the n-correlated analysis we briefly turn to the quantities ^6,7 defined in 
eqs. ffity and A plateau in the corresponding effective-C6,7 plots for a particular pair 
{n,rn} at a value inconsistent with zero would signal the breakdown of the large-volume 
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FIG. 6: Effective m^ii 2) plots for the LO (top-left), NLO (top-right) and N 2 LO (bottom-left) 
extractions using the n = 8 energy and the N 3 LO effective m^a^ 2 ^ extracted from {n, m} = {8, 4} 



(bottom- right) from eq. (10). In each case the m n = 291 MeV ensemble is used and the scale in 
each plot is identical. The horizontal band corresponds to the extraction using the exact eigenvalue 
method for n = 2. 



expansion in eq. ([T]) that is central to our analysis. In all cases, no such breakdown is seen. 
However for increasing n and m, the uncertainties increase. For example, for the systems 
with {n, m} = {11, 12} at = 291 MeV, these quantities are found to be C6 = C7 = 0.0(3). 



E. n-correlated analysis 



The most complete use of the full set of energy differences that we have computed is made by 
performing the coupled, 0(L~ 7 ) analysis of the n — 2, ... ,12 effective energy differences to 

extract a^ 2 ^ and Tj 3 , including the correlations in both t and n as discussed in the preceding 



section. The resultant fits of such an analysis are shown in figs. [TT]-[T4]for the four ensembles. 
The extracted fit parameters, a^ 2 ^ and rj 3 , are central results of this work. They are given 

For comparison with the 



in Table III and their uncertainty ellipses are shown in fig. [15| 

simpler analysis above, the shaded regions in fig. [lO correspond to the values of rfe extracted 
from this correlated analysis and are seen to be consistent with the extractions made using 
eq. (11). Similarly, the extracted a^ 2 -* agrees with that obtained from the exact eigenvalue 
2 (and hence with all the N 3 LO extractions above). 



method for n 
In fig. 16 



7/3 is plotted versus m^j /„., in units of the estimate based upon NDA, %' Nj 1 



l/(m,r/^), as discussed in Ref. [5], [8] . While the three-body interaction is consistent with its 
NDA estimate for the lightest three pion masses, it is found to be consistent with zero at the 
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_(J=2) 

FIG. 7: Extractions of m^a^ for each ensemble. The blue circles, mauve squares and brown 
diamonds correspond to the extractions at 0(L -3 ), C(L -4 ), and 0(L~ 5 ) in the 1/L expansion 
given in eq. Q. The stars correspond to the 0(L -6 ) extractions using eq. (10), which requires two 
different energy shifts. At any given n, we have shown various combinations of m < n. Finally the 

_(I=2) 

solid line and the shaded region correspond to the extraction of m^a\^ using the exact eigenvalue 
method of Liischer from the n = 2 data. In all cases, statistical and systematic uncertainties 
have been combined in quadrature. The upper-left, upper-right, lower-left and lower-right panels 
corresponds to m n = 291,352,491,591 MeV, respectively. 



TABLE III: The results of simultaneous fit of and ry 3 to all twelve correlation functions, 

( 1=2) 

n = 2,3, ..,12, on each ensemble. (Recall that the nuclear physics sign convention for a^^ is 
being used.) 



Quantity 


~ 291 MeV 


m* ~ 352 MeV 


Hi, ~ 491 MeV 


77V ~ 591 MeV 




1.990(11)(14) 


2.3230(57) (30) 


3.0585(49) (95) 


3.4758(98) (60) 


_(J=2) 


0.1644(40) ( _ + « ) 


0.2058(45) ( +2 ) 


0.3497(69) ( ) 


0.4761(96) ( + ™) 




0.80(09) ( ) 


1.02(08) ( 


0.90(12) 


0.55(23) (™) 


X 2 /dof (dof) 


1.3 (65) 


1.9 (63) 


1.8 (63) 


1.3 (53) 



heaviest mass, = 591 MeV. It is desirable to reduce the uncertainties in this calculation 
to see, if in fact, the three-body interaction is decreasing with increasing pion mass. If this is 
found to be the case then a more detailed study in this high pion mass region is warranted. 

It is interesting to study how the larger n energy differences influence the extraction of 
the two parameters. To do so, we have performed a series of fits including only the energy 
differences up to a given n max . The resulting confidence regions for the parameters extracted 
from the m n = 291 MeV ensemble are shown in fig. [l7]for n max = 3, 5, 7, 9, 11, 12. Clearly 
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FIG. 8: The effective energy splitting plots for the 



291 MeV ensemble. The solid lines 
l^ 2 ^ from the tt + it^ 



energy- 



correspond to the energy differences of eq. using the value of a 
splitting (using the exact eigenvalue method) and setting 773 = 0. The shaded region shows the 
statistical and systematic uncertainties combined in quadrature. 



the inclusion of higher n data improves the determination of r] 3 in particular. 
V. THE THREE-PION INTERACTION, ^ , rj%, AND r? 3 (/i) 

In the preceding Section the 7r + 7r + 7r + -interaction, 773 , has been extracted successfully from 
the lattice QCD calculations and is non-zero at the lighter quark masses used in this study. 
The renormalization group invariant, but volume dependent three-body interaction that 
arises at 0(L~ 6 ), rj%, that receives perturbative corrections in the volume expansion to 
become % at 0(L~ 7 ), is somewhat more problematic to extract from these particular lattice 
calculations due to the size of a^ 2 ^ in relation to the size of the lattice. 

There are two terms that must be calculated in order to determine r\\ , one is additive 
and one is multiplicative, as given in eq. (|4]). Assuming that the effective range for 7r + 7r + 
scattering is not orders of magnitude larger than afJ= 2 \ the additive term makes a very 
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FIG. 9: The effective r/ 3 plots for the = 291 MeV ensemble using eq. rtllfc at 0{L- 7 ). Statistical 
and systematic uncertainties are added in quadrature (shaded band). 



small contribution to r] 3 . To make this explicit, it is useful to rewrite eq. (4) as 

m * ft % = a m m * ft V3 + A» = » 

where 



a. 



'/■-! 



(24) 



(25) 



The numerical values of (3 m are shown in Table IV, and are all seen to be very small for 
r ~ \a\. The multiplicative factor, a V3 , is also shown in Table IV, and its values make 



clear that the perturbative expansion for the three-body interaction is converging extremely 
slowly at the lightest pion mass, with a ~ 75% correction at 0(L~ 7 ) to the 0(L~ 6 ) result, 
and fails completely at both = 491 MeV and 591 MeV for the lattice volumes used in 
this work. Clearly larger lattice volumes will be required in order to determine the three- 
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FIG. 10: The value extracted for r/3 using eq. (jllj) at 0(L 7 ) as a function of n for each ensemble 
used in the calculation. The horizontal line in each plot corresponds to the value extracted using 



all the data in the n-correlated fit of Section IV E The inner and outer shaded bands correspond to 



the statistical uncertainty and the statistical and systematic uncertainties combined in quadrature, 
respectively. 

body interaction 7J3 with precision. 14 This is in stark contrast to the two-body scattering 
parameters which can be extracted to high precision from these same lattice volumes and 
the perturbative expansion in eq. ([I]) appears to be converging. So while it is possible to 



=L 



TABLE IV: Correction factors required to determine 773 from r/ 3 , defined in eq. (25). 



Quantity 


m, ~ 291 MeV 


in, ~ 352 MeV 


in, ~ 491 MeV 


m w ~ 591 MeV 




1.990(11)(14) 


2.3230(57)(30) 


3.0585(49) (95) 


3.4758(98) (60) 




1.74(3) 


1.78(2) 


1.97(2) 


2.08(3) 


An 


-0.0038(7) 


-0.0056(7) 


-0.020(2) 


-0.044(6) 


a„ 3 [(3.5 fm) 3 ] 


1.53(2) 


1.56(1) 


1.69(2) 


1.77(2) 


A, 3 [(3.5 fm) 3 ] 


-0.0027(5) 


-0.0040(5) 


-0.014(2) 


-0.031(4) 



determine the "dressed" three-pion interaction, the uncertainty in the determination of the 
bare three-pion interaction is large, not because of the uncertainty in the lattice calculation, 
but due to the relatively large higher-order terms in the volume expansion when evaluated at 



this present lattice volume (a theoretical systematic uncertainty). In Table IV, the correction 
factors are also given for a lattice of spatial volume (3.5 fm) 3 for comparison. However, 
even in these large volumes the correction factors at one higher order in the large volume 



We note that it is t] 3 that was extracted in Ref. [3]. 
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FIG. 11: n-correlated fits to the n = 2, . . . , 12 energy differences for the = 291 MeV ensem- 
ble. All energy differences, AE^g(t) for n = 2, ... ,12, are used. Statistical, and statistical plus 
systematic (added in quadrature) uncertainties are shown as the inner and outer shaded regions, 
respectively. 



expansion of ~ 50% and clearly even larger volumes, L > 3.5 fm, are required in order to 
have the volume expansion of the three-pion interaction under perturbative control. 

The quantity 773 ([/,) is a renormalization scheme dependent quantity that is independent 
of the volume, and as such is the quantity that most directly enters into the calculation of 
other many-body processes. It is easily extracted from r\\ via eq. (IHJ. However, given the 
present theoretical systematic uncertainties in 7/3, we do not attempt a determination of 



VI. THE EQUATION OF STATE AND THE ISOSPIN CHEMICAL POTENTIAL 

The energy of the n-ix + system as a function of volume and of the number of 7r + 's is given 
explicitly in eq. (TO) in the large-volume expansion. From the equation of state, the isospin 
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FIG. 12: ra-correlated fits to the n = 2, . . . , 12 energy differences for the m, = 352 MeV ensemble. 



chemical potential is defined as 

Hi 



dE n 



dn 



(26) 

y=const 



which can be constructed analytical from eq. Q or numerically from the results of the 
lattice calculation by using a simple finite difference approximation. The resulting ratio 
of this isospin-chemical potential to the pion mass for each of our ensembles is shown in 



fig. 18 as a function of n, and for convenience, the isospin density, pj (which in this case 
is the number density). We note that for the 12-7r + system, the number density is p\ 



12/ L 3 = 0.77 far 3 . In fig. Il8l the solid curves corresponds to the prediction at 0(L 7 ) from 
eq. (jl]), and this prediction differs insignificantly from that at 0(L~ 6 ). There have been 
recent works that perform lattice QCD calculations at finite isospin chemical potential, e.g. 
Ref. [62, 63J, where the starting point is a chemical potential term for the quarks added to 
the QCD action. Compelling results for the formation of the charged-pion condensate at 
Hi ~ m-x have been produced. 

Pionic systems at finite isospin chemical potential (and temperature) have been explored 
theoretically in x?T j6U [651 ES |67], primarily as a step toward understanding finite den- 
sity nuclear systems. The inclusion of the isospin chemical potential can be accomplished 
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FIG. 13: re-correlated fits to the n = 2, . . . , 12 energy differences for the m n = 491 MeV ensemble. 



straightforwardly by replacing the time- component of the covariant derivative acting on the 
exponential field of 7r's, DqT,, with VoS = DqT, — ipi\ [ra , £ ]. For /// < the vacuum 
state is the same as it is for pj = 0, S = 1, but for //j > the vacuum alignment changes 
and becomes, 

S = cos a + i (ti cos + r 2 sin 0) sin a , (27) 

using the notation of Refs. [661 EZ]- Minimizing the potential energy at LO in gives 
cos a = m^j p\ and a relation between the isospin density, pj, and chemical potential 15 , 

Pi = f 1-^1 = 2/. 2 (^-m.) - 3^ (^-m.) 2 + ... . (28) 



\ pj J m 

By construction, our lattice QCD calculation is in the regime of pi > m n , and we expect 
that these LO results should come perturbatively close to describing the properties of 



15 



The numerical factors that appear in eq. (28 1 differ by factors of two from Refs. [64] [65j [66l |67j due to 



the definition of /„■. At the physical pion mass we use f„ ~ 132 MeV. 
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FIG. 14: ra-correlated fits to the n = 2, . . . , 12 energy differences for the m n = 591 MeV ensemble. 



the n-7t + systems. Implicit in this LO xPT analysis are not only the 7r + 7r + -interactions, but 
also multi-pion interactions, including the 7r + 7r + 7r + -interaction. In fact, the NDA estimate 
of the 7r + 7r + 7r + -interaction arises from the LO terms in %PT. In fig. 19, we show the isospin 
chemical potential versus the number of 7r + 's in the lattice volume, which can be directly 
translated into isospin density, pj = n/L 3 . To determine the importance of the 7r + 7r + 7r + - 
interaction, we show the \ii calculated directly from the lattice calculation. We also show 



the curve resulting from eq. ( l| with the values of 2 -* and r/ 3 extracted from the lattice 
calculation and their associated correlated uncertainties (red curves and shaded regions), 
and we show the curve resulting from setting rj%= in eq. |IJ (blue curve and shaded 
region). It is clear from fig. 19 that the 7r + 7r + 7r + -interaction plays an important role in the 



relation between the isospin chemical potential and the isospin density. Further, in fig. 19 



the dashed curve is the result of LO x?T, as given in eq. (28), which is seen to describe the 
result of the lattice calculation well at all the pion masses we have explored. There does 
appear to be a slight m^-dependent systematic difference, but the magnitude of this effect 
is consistent with terms higher order in x?T. 

Another quantity of interest that one wishes to determine for such systems is the pressure 
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FIG. 15: The 68% and 90% confidence regions for the two parameter fits performed to the n = 
2, ... ,12 energy differences for each ensemble used in our work. Only statistical uncertainties are 
shown. 



as a function of the pj, 



1 dE r . 



3L 2 dL 



(29) 



n=const 



Unfortunately, we cannot recover the pressure directly from our lattice calculations using 
eq. (jl) because both a^ 2 ^ and rj^ depend implicitly upon the volume. Further, without 
lattice calculations of different volumes, the pressure cannot be determined directly from 
the lattice calculations either. Therefore, with the present lattice calculations we have 
performed, while the isospin chemical potential can be determined as a function of density, 
the pressure cannot. 

An important issue to consider is the impact of the finite lattice volume upon the relation 
between the isospin chemical potential and the isospin density. The periodic boundary 
conditions imposed in the spatial directions leave the zero-mode untouched, but discretize 
the non-zero modes. As such, the tree-level matrix element of the interaction Hamiltonian 
between initial and final states in which each pion carries zero three-momentum is unaffected 
by the finite volume. However, the boundary conditions will modify contributions at the 
one-loop level and beyond. Therefore, in the limit where the contributions from loop- 
level diagrams are small, as is the case for a small scattering length, the modifications 
of the relation between isospin chemical potential and isospin density due to the boundary 
conditions is expected to be perturbatively small 16 . In systems with large scattering lengths, 



16 



In infinite volume, the leading contributions to the ground state energy of an imperfect dilute Bose gas 
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e.g. nucleons, the finite-volume modifications may be substantial, and this requires further 
study. 



are calculated easily by performing a Bogoliubov transformation on the Hamiltonian, leaving the leading 
order mean-field contribution and a term resulting from a summation over the non-zero momentum states. 
The leading terms in the density expansion of the energy per particle on the ground-state are 



E/N 



2irap 
M 



(30) 



When placed in a finite cubic volume with periodic boundary conditions, the integral over non-zero 
momentum states that gives rise to the term that is non-analytic in the number-density, p, is replaced by 
a summation over the allowed momentum states in the volume. When the unperturbed level spacing in 
the volume is large compared with the energy-shift at leading order in the density expansion, the contents 



of the square brackets in eq. ( 30 1 becomes 




1 - 


(A) 1 + ■■■] 







(31) 



recovering the leading terms in the finite- volume expansion in eq. ([I]) as the number of particles becomes 
large. See also Ref. [9]. 
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FIG. 17: Dependence of the fit parameters, and f] 3 , on the number of energy differences 

included in the fit for the m n = 291 MeV ensemble, n < n max implies that the energies up to n max 
are included in the correlated fit. The inner and outer ellipses show the 68% and 90% confidence 
regions. 



VII. CONCLUSIONS 

One of the major challenges facing nuclear physics is the solution of the hadronic many-body 
problem, including the effects of the multi-hadron interactions induced at the scale of chiral 
symmetry breaking. While the two-nucleon interaction is overwhelmingly the dominant 
interaction in multi-nucleon systems, a three-nucleon interaction contributes significantly to 
the structure and interactions of nuclei. Calculating the properties and interactions of nuclei 
is a major goal of lattice QCD and as a small step in this direction we have performed the 
first lattice QCD calculations of the properties of systems comprised of multiple-7r + 's. The 
present paper is a detailed follow-up to the letter we recently published [I], and extends 
that work from studies of systems comprised of up to five pions, to systems containing up 
to twelve pions. As a technical detail, this required calculations of very high (64 decimal 
digit) precision. 

We have convincingly determined a non-zero value for the dressed three-body interaction, 
rjg, defined in eq. (fll), at the lightest three pion masses, = 291,352 and 491 MeV, and 
find a value consistent with zero at the heaviest mass, = 591 MeV. These central 



results are summarized in Table III and fig 16 above. Given the lattice volumes in which 



the calculations have been performed, and the value of the 7r + 7r + scattering length, the 



connection between rj^ and the underlying three-body interaction, 773 (/x), (see eqs. (4) and 
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FIG. 18: The isospin chemical potential as a function of the number of pions in a fixed volume 
(equivalent to the isospin density pj) . The points and their associated uncertainties are the results 
of the lattice calculation where a finite-difference has been used to construct the derivative with 
respect to the number of 7r + 's. The solid curves and bands result from the analytic expression for 
the energy of the ground state in the large volume expansion, eq. using the fit values for a^~ 2 ^ 
and 773 an d their correlated uncertainties. 

(j5|) is slowly converging, making a meaningful extraction of 773 (/x) currently impractical. 
Clearly, in addition to higher order calculations (0(L~ 8 ) and beyond) of the energy of 
multi-pion systems in a finite volume, lattice calculations in larger volumes, and of excited 
states in these volumes, are required to disentangle the bare from the dressed interactions. 

An important outcome of our calculations is a determination of the relation between 
the isospin chemical potential and the isospin density. For m n < 400 MeV the 7r + 7r + 7r + 
interaction is seen to make a sizable contribution to this relation. Such contributions are 
implicit in the LO %PT theoretical result, but our calculation has provided the first explicit 
QCD calculation of the importance of multi-pion interactions. 

This work and our previous Letter represent the first study of multi-hadrons systems di- 
rectly from QCD. Whilst encouraging, we conclude by reiterating that a significant amount 
of theoretical work remains to be performed in order to explore meaningfully more compli- 
cated systems such as those involving large numbers of baryons. 
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the energy of the ground state in the large volume expansion, eq. (lb, using the fit values for a 
j — 

and r/ 3 and their correlated uncertainties. The solid blue (darker) curves and bands are similarly 
the results for the fitted 
order prediction of xPT. 
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APPENDIX A: CONTRACTIONS FOR I z = n < 13 MESONS 

The contractions required for the n — 1,. . . , 13 correlation functions are given below. The 
notation 

Tj = tr c ,s [IF] 

is used for brevity with the trace being over color and spinor indices. The index, j, is the 



matrix power to which II is raised. 

= Ti, (Al) 

C 2 (t) = I?-T 2) (A2) 

C 3 (t) = T? - 3T 2 Ti + 2T 3 , (A3) 

C 4 (i) = 7? - GTaTf + 8T 3 Ti + 3T| - 6T 4 , (A4) 

Cs(i) = -\QT 2 Tf + 2QT 3 Tl + lbT 2 2 T l -?,QT 4 T 1 -2QT 2 T 3 + 2AT bl (A5) 

C 6 (t) = rf - 15T 2 T? + 40T 3 T? + 45T 2 2 T 1 2 - 90T 4 T? - 120T 2 T 3 Ti 

+144T 5 Ti - 15T 2 3 + 40T 3 2 + 90T 2 T 4 - 120T 6 , (A6) 

C 7 (t) = Tl - 2\T 2 T\ + 70T 3 T? + lOST^T? - 210T A T? - 420T 2 T 3 T? + 504T 5 T? - 105T 2 3 Ti 

+280T|Ti + 630T 2 T 4 Ti - 840T 6 Ti + 210T 2 2 T 3 - 420T 3 T 4 - 504T 2 T 5 + 720T 7 , (A7) 

C 8 (*) = Tf - 28^^ + mi^Tf + 210T 2 2 7 1 4 - 420T 4 T 4 - 1120T 2 T 3 T? + 1344T51; 3 - 420T 2 3 T 1 2 

+ 1120T 3 2 T 1 2 + 2b2QT 2 T A Tl - 3360T 6 T? + 1680T 2 2 T 3 Ti - 3360T 3 T 4 Ti - 4032T 2 T 5 Ti + 5760T 7 Ti 

+ 105T 2 4 - 1120T 2 r 3 2 + 1260T| - 1260T 2 2 T 4 + 2688T 3 T 5 + 3360T 2 T 6 - 5040T 8 , (A8) 



C 9 (t) = Tl - 36T 2 T7 + 168^^ + 378T 2 2 T 1 5 - 756T 4 T-f - 2520T 2 T 3 T 1 4 + 3024T 5 T 4 - 1260T 2 3 Tf 
+3360T 3 2 T 1 3 + 7560T 2 T 4 T 1 3 - 10080^^ + 7560T 2 2 T 3 7 1 2 - 15120T 3 T 4 T? - 18U4T 2 T 5 T? 
+25920T 7 T? + 945T 2 4 Ti - 10080T 2 T 3 2 Ti + 11340T 4 2 T! - 11340T 2 2 T 4 Ti + 24192T 3 T 5 Ti 
+30240T 2 T 6 Ti - 45360T 8 Ti + 2240T 3 3 - 2520T 2 3 T 3 + 15120T 2 T 3 T 4 + 9072T 2 2 T 5 - 18144T 4 T 5 
-20160T 3 T 6 - 25920T 2 T 7 + 40320T 9 , (A9) 

C10O) = T i° - 45X 2 T 1 8 + 2AQT 3 Tl + 630T 2 2 7f - 1260T 4 T-f - 5040T 2 T 3 T-f + 6048T 5 T? - 3150T 2 3 T 4 
+8400T|T 4 + 18900T 2 T 4 T 4 - 25200^^ + 25200T 2 2 T 3 T 1 3 - 50400T 3 T 4 T? - 60480T 2 r 5 ri 5 
+86400T 7 T 3 + 4725T 2 4 T 2 - 50400T 2 T 2 T 2 + 56700T 4 2 T 2 - 56700T 2 2 T 4 T 2 + 120960T 3 T 5 T 2 
+151200T 2 r 6 r 2 - 226800^^ + 22400T 3 3 Ti - 25200T 2 3 T 3 Ti + 151200T 2 T 3 T 4 Ti + 90720T 2 2 T 5 Ti 
-181440T 4 T 5 Ti - 201600T 3 r 6 Ti - 259200T 2 T 7 Ti + 403200T 9 Ti - 945T 2 5 + 25200T 2 2 T 3 2 
-56700T 2 T 2 + 72576T 2 + 18900T 3 T 4 - 50400T 2 T 4 - 120960T 2 T 3 T 5 - 75600T 2 2 T 6 
+151200T 4 T 6 + 172800T 3 T 7 + 226800T 2 T 8 - 362880Ti , (A10) 
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Cu(t) = T 11 - 55T 2 T? + 330T 3 T? + 990T$T? - 1980T 4 T^ ~ 9240T 2 T 3 T? + llOSSTgTf - 6930T$T? 

+18480T 3 2 T 1 5 + 41580T 2 r 4 ri i - 55440^^ + 69300T 2 2 T 3 7^ - 138600T3T4T; 4 - 166320T 2 T 5 T^ 
+237600T 7 T 4 + 17325T 2 4 T 1 3 - 184800T 2 T 3 2 T 3 + 207900T 4 2 T 3 - 207900T 2 2 T 4 T 3 + 443520T3T5T 3 
+554400^^^ - 831600787? + 123200T 3 3 T 1 2 - 138600T$ T 3 T? + 83imQT 2 T 3 T A Tl 
+498960T 2 2 T 5 T 1 2 - 997920T4T5T! 2 - llOSSOOTgTeT! 2 - 1425600T 2 T 7 T? + 2217600T9T; 2 
-10395T 2 5 Ti + 277200T 2 2 T 3 2 Ti - 6237007 2 T 4 2 Ti + 798336T 5 2 7i + 207900T 2 3 T 4 Ti - 5544007 1 3 2 T 4 7i 
-1330560r 2 T 3 T 5 Ti - 831600T 2 2 T 6 Ti + 1663200T 4 T 6 7i + 1900800T 3 T 7 Ti + 2494800T 2 T 8 Ti 
-3991680Ti Ti - 123200T 2 T 3 3 + 415800T 3 T 4 2 + 34650T 2 4 T 3 - 415800T 2 2 T 3 T 4 - 166320T 2 3 T 5 
+443520T 3 2 r 5 + 997920T 2 T 4 T 5 + 1108800T 2 T 3 T 6 - 1330560T 5 T 6 + 712800T 2 2 T 7 - 1425600T 4 T 7 
-1663200T 3 T 8 - 2217600T 2 T 9 + 3628800Tu , (All) 

Ci 2 {t) = 7\ 12 - 66T 2 T 10 + 440737? + 1485717? - 2970T 4 T? - 15840^^^ + I9OO8T5I? - 13860T 2 3 7? 

+369607f7? + 83160T 2 T 4 T? - 110880^^ + 1663207 2 2 T 3 T 1 5 - 332640T3T47? - 399168T 2 TgT? 
+570240T 7 I? + 51975T.f7? - 5544007 1 2 T 2 T 4 + 623700T 4 2 T 4 - 623700T 2 2 T 4 T 4 + 1330560T 3 r 5 r 4 
+1663200^^^ - 2494800T8T] 4 + 492800T|7? - 554400T 2 3 T 3 T 1 3 + 3326400T 2 T 3 T 4 7? 
+1995840T 2 2 T 5 T 1 3 - 3991680T 4 T57? - 4435200T 3 T 6 7? - 5702400T 2 T 7 7? + 887040071,7? 
-62370T 2 5 T 2 + 1663200T 2 2 T 2 T 2 - 37422007 2 T 4 2 T 2 + 4790016T 2 T 2 + 12474007 , 2 3 7 4 T? 
-3326400T 3 2 T 4 T 1 2 - 7983360T 2 T 3 7 1 5 7? - 4989600T 2 2 T 6 T 1 2 + 9979200T4T67? + 11404800T 3 r 7 T 1 2 
+14968800^^^ - 239500807i 7? - 1478400T 2 T 3 3 7i + 498960073T 4 2 7i + 415800T 2 4 T 3 Ti 
-4989600T 2 2 T 3 T 4 Ti - 1995840T 2 3 T 5 7i + 5322240T 3 2 T 5 7i + 11975040T 2 7 4 7 1 5 7i + 13305600T 2 T 3 T 6 Ti 
-15966720T 5 T 6 Ti + 8553600T 2 2 T 7 7i - 17107200T 4 T 7 Ti - 19958400T 3 T 8 Ti - 26611200T 2 T 9 Ti 
+43545600T n Ti + 10395T 2 6 + 246400T 4 - 12474007? - 5544007?T 3 2 + 1871100T 2 2 7? 
-4790016T 2 7f + 6652800T 6 2 - 311850T 2 4 T 4 + 33264007 2 T 3 2 7 4 + 3991680T 2 2 r 3 T 5 - 7983360T 3 T 4 T 5 
+16632007?T 6 - 44352007 1 3 2 T 6 - 9979200T 2 T 4 T 6 - 1 1404800T 2 T 3 T 7 + 13685760T 5 T 7 
-7484400T 2 2 T 8 + 14968800T 4 T 8 + 17740800T 3 T 9 + 23950080T 2 T 10 - 399168007\ 2 , (A12) 
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C 13 (t) = T 13 - 78T2TI 1 + 572T 3 T 10 + 2145T 2 2 T 1 9 - 4290T 4 T? - 257iOT 2 T 3 T? + 30888T s Tf 
-25740lf 7? + 686407f T- 7 + 154440T 2 T 4 T 7 - 205920T 6 T 7 + 360360T2t 3 T? 
-720720T 3 TiT^ - 864864T 2 Tglf + 1235520T7I; 6 + 135135T 2 4 T 1 5 - l4A144QT 2 T£lf 
+1621620T 4 2 T 1 5 - 1621620T 2 2 T^Tf + 3459456T3T57; 5 + 4324320T 2 T 6 Tf - 64864807 8 7 T 1 5 
+1601600T 3 3 T 1 4 - 1801800Tf TgT- 4 + IO8IO8OOT2T3T4T 4 + 6486480T 2 2 T 5 T 4 - 12972960T 4 !T 5 T 4 
-144144007 3 7' 6 7 1 4 - 18532800T 2 T 7 T^ + 28828800T 9 T 4 - 270270T 2 5 T 1 3 + 7207200T 2 2 T 3 2 T 1 3 
-16216200T 2 T|l? + 20756736T 5 2 T 1 3 + 5405400T 2 3 T 4 T 1 3 - 14414400T 3 2 T 4 r l 3 - 34594560T 2 TaTsi; 3 
-2162160071^^ + 43243200T 4 T 6 !T 1 3 + 4942080073^^ + 64864800T 2 7 8 T 1 3 - 103783680^0^ 
-9609600T 2 T$T? + 32432400T 3 T 4 2 T 1 2 + 2702700T 2 4 T 3 T 1 2 - 32432400T 2 2 T 3 T 4 7 1 2 
-12972960T 2 3 T 5 T 1 2 + 34594560T 3 2 T 5 T 1 2 + 77837760T 2 T 4 T57\ 2 + 86486400T 2 T 3 T 6 T 1 2 
-103783680T5T6T! 2 + 55598400T 2 2 T 7 T 1 2 - IIII968OOT4T77? - 129729600T 3 T 8 Tf 
-172972800T 2 T 9 T 1 2 + 283046400^1^ + 135135T 2 6 Ti + 3203200T 3 4 Ti - 16216200r 4 3 T! 
-7207200T 3 T 2 T 1 + 24324300T 2 2 T 4 2 T 1 - 62270208T 2 T 2 T 1 + 86486400T 2 T 1 
-4054050T 2 4 r 4 Ti + 43243200T 2 T 3 2 T 4 Ti + 51891840T 2 2 T 3 r 5 ri - 103783680T3T 4 T 5 Ti 
+21621600T 2 3 T 6 Ti - 57657600T 3 2 T 6 7i - 12972960Q7 , 2 7 4 T 6 7 1 - 148262400T 2 T 3 T 7 Ti 
+177914880T5T7T! - 97297200T 2 2 T 8 Ti + 194594400^^^ + 230630400T 3 T 9 Ti 
+311351040T 2 7 1 io7 n i - 518918400Ti 2 Ti + 4804800T 2 2 T 3 3 - 32432400T 2 T 3 T 4 2 
+41513472T 3 r 5 2 - 540540T 2 5 T 3 - 9609600T 3 3 T 4 + 10810800T 2 3 T 3 T 4 
+3243240T 4 T 5 - 34594560T 2 T 3 2 T 5 + 38918880T 2 T 5 - 38918880T 2 2 T 4 T 5 
-43243200T 2 2 r 3 T 6 + 86486400T 3 T 4 T 6 + 103783680T 2 T 5 T 6 - 18532800T 2 3 T 7 
+49420800T 3 2 T 7 + 111196800T 2 T 4 T 7 - 148262400T 6 T 7 + 129729600T 2 T 3 T 8 
-155675520r 5 T 8 + 86486400T 2 2 T 9 - 172972800T 4 T 9 - 207567360T 3 Ti 

-283046400r 2 Tn + 4790016007i 3 . (A13) 



APPENDIX B: NUMERICAL PRECISION IN n-MESON CORRELATION 
FUNCTIONS 

In order to extract the correlation functions for systems with n > 8 mesons it is necessary 
to perform the contractions using high precision numerical techniques. It is necessary to 
calculate the propagators to an analogous level of precision. 



1. Contraction Precision 



As discussed in the main text, the double precision numerical representations are insufficient 
to accurately compute the propagator contractions required for the large n correlators. Here 
we exhibit the failure of these calculations and discuss their origin. 



In figs. [20] and |2TJ we 

compare the n particle effective energies and effective energy differences performed using 
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FIG. 20: Comparison of double precision contractions and 64 decimal digit precision contractions 
for the effective energy plots for n = 9, 10, 11 and 12 pions. 
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FIG. 21: Comparison of double precision contractions and 64 decimal digit precision contractions 
for the effective energy difference plots for n = 9, 10, 11 and 12 pions. 



double precision (squares) and using 64 decimal digit internal precision (stars, slightly offset 
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to the right for clarity) for n = 9, 10, 11 and 12 for the = 352 MeV ensemble. The 
breakdown of double-precision calculations with increasing n and t, where non-zero values 
of the correlation functions are lost to noise, is clear. This breakdown has its origin in 
the large powers to which propagators are raised in the contractions required for systems 
with a large number of mesons. Viewing the zero three-momentum propagator as a time- 
dependent 12 x 12 matrix, the range of numbers that can be represented in double precision 
is too limited and small elements of high powers of these propagators are rounded away 
in tracing or other matrix manipulations. Such terms are crucial for maintaining gauge 
invariance, and their removal leads to a serious degradation of the correlation function. For 
the numbers of mesons considered in this work, n < 13, it is sufficient to use 64 decimal 
digit precision, but for n > 13, such precision will rapidly become insufficient. 



2. Propagator Precision 



Given a particular gauge field, the accuracy with which the quark propagator is computed 
(e.g. the tolerance in conjugate gradient (CG) algorithm) will influence the correlation 
functions computed from that propagator. In extreme cases, where the solution is quite 
inaccurate in comparison to the statistical precision of the importance sampling procedure, 
this residual error can persist through the ensemble averaging and lead to errors in the 
correlation functions. In the case of the multi-pion correlations studied here, the situation 
is particularly acute as the high powers to which the propagators are taken can enhance 
small effects. At this point in time we have not performed a detailed study of this issue. It 
would require generating full sets of propagators at different inversion precisions on the same 
configurations. We have simply studied the effect on a few representative configurations, 
and it is clear that for large n, significant differences in the correlation functions can arise 
from loose tolerance of the CG-solver. In future studies of even larger n (and very high 
precision studies of any observable), this issue must be investigated in more detail. 
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